function ShiftedLognormal
clc;close all;

% Monte Carlo settings

NoOfPaths  = 5000;
NoOfSteps  = 1000; 

% Libor parameters

T         = 3.0;
sigma     = 0.2;
L0        = -0.05;
shift     = 0.1;
K         = [0.95];
CP        = 'c';
P0T       =@(T)exp(-0.1*T);

% Shifted path generation

[L,time] = GeneratePathsGBMShifted(NoOfPaths,NoOfSteps,T,0.0,sigma,L0,shift);
figure(1)
plot(time,L(1:25,:))
grid on
xlim([0,T])
xlabel('time')
ylabel('Libor L(t)')
title('Shifted Lognormal Paths')

% Shifted lognormal density

f_L = @(x,shiftArg,T)lognpdf(x+shiftArg,log(L0+shiftArg)+(-0.5*sigma^2)*T,sigma*sqrt(T));
shiftV = [1.0, 2.0, 3.0, 4.0];
argLegend = cell(length(shiftV),1);
nGrid = 50;
xArgM=zeros(length(shiftV),nGrid);
pdfM =zeros(length(shiftV),nGrid);
idx = 1;
for shiftTemp=shiftV
    x =linspace(-shiftTemp,6,nGrid);
    f_shifted = f_L(x,shiftTemp,T);
    xArgM(idx,:)=x;
    pdfM(idx,:)=f_shifted;
    argLegend{idx} = sprintf('shift = %.1f',shiftTemp);
    idx = idx +1;
end
MakeFigure(x, pdfM, argLegend,'','shifted density')

% Paths and 3D graph

figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');

% Density plot for S(t)

plot(time,L(1:25,:),'linewidth',1,'color',[0 0.45 0.74])
hold on
TVec=linspace(0,T,5);
for i=1:length(TVec)
    ind=find(time>=TVec(i),1,'first');
    x_arg=linspace(min(L(:,ind))-L0-shift,max(L(:,ind))+L0+shift,100);
    plot3(TVec(i)*ones(length(x_arg),1),x_arg,f_L(x_arg,shift,TVec(i)),'k','linewidth',2)
end
axis([0,TVec(end),-shift,max(max(L(1:25,:)))])
grid on;
xlabel('t')
ylabel('$$\ell(t)$$')
zlabel('density')
view(axes1,[-68.8 40.08]);

% Caplet/Floorlet: Monte Carlo vs. analytic expression

K = linspace(-shift,abs(L0)*3,25);
optPrice = zeros(length(K),1);
idx = 1;
for k=K
    optPrice(idx)= P0T(T)* EUOptionPriceFromMCPaths(CP,L(:,end),k,T,0.0);
    idx = idx +1;
end

% Analytic expression

optPriceAnalytical = P0T(T)*BS_Call_Put_Option_Price_Shifted(CP,L0,K,sigma,T,0.0,shift);
figure(4)
plot(K,optPrice,'-b','linewidth',2.0)
hold on
plot(K,optPriceAnalytical,'--r','linewidth',2.0)
grid on
legend('Option Price: Monte Carlo','Option Price: Analytical')


function optValue = EUOptionPriceFromMCPaths(CP,S,K,T,r)
if lower(CP) == 'c' || lower(CP) == 1
    optValue = exp(-r*T)*mean(max(S-K,0));
elseif lower(CP) == 'p' || lower(CP) == -1
    optValue = exp(-r*T)*mean(max(K-S,0));
end

function [S_shifted,time] = GeneratePathsGBMShifted(NoOfPaths,NoOfSteps,T,r,sigma,S_0,shift)
S0_shifted = S_0 + shift;
[S,time] = GeneratePathsGBMEuler(NoOfPaths,NoOfSteps,T,r,sigma,S0_shifted);
S_shifted = S - shift;

function [S,time] = GeneratePathsGBMEuler(NoOfPaths,NoOfSteps,T,r,sigma,S0)

% Approximation

S=zeros(NoOfPaths,NoOfSteps);
S(:,1) = S0;

% Random noise

Z=random('normal',0,1,[NoOfPaths,NoOfSteps]);
W=zeros([NoOfPaths,NoOfSteps]);

dt = T / NoOfSteps;
time = zeros([NoOfSteps+1,1]);
for i=1:NoOfSteps
    if NoOfPaths>1
        Z(:,i)   = (Z(:,i) - mean(Z(:,i))) / std(Z(:,i));
    end
    W(:,i+1)  = W(:,i) + sqrt(dt).*Z(:,i);
    S(:,i+1) = S(:,i) + r * S(:,i)*dt + sigma * S(:,i).*(W(:,i+1)-W(:,i));
    time(i+1) = time(i) + dt;
end

function value = BS_Call_Put_Option_Price_Shifted(CP,S_0,K,sigma,tau,r,shift)
K_new = K + shift;
S_0_new = S_0 + shift;
value = BS_Call_Put_Option_Price(CP,S_0_new,K_new,sigma,tau,r);

% Closed-form expression of European call/put option with Black-Scholes formula

function value=BS_Call_Put_Option_Price(CP,S_0,K,sigma,tau,r)
d1    = (log(S_0 ./ K) + (r + 0.5 * sigma^2) * tau) / (sigma * sqrt(tau));
d2    = d1 - sigma * sqrt(tau);
if lower(CP) == 'c' || lower(CP) == 1
    value =normcdf(d1) * S_0 - normcdf(d2) .* K * exp(-r * tau);
elseif lower(CP) == 'p' || lower(CP) == -1
    value =normcdf(-d2) .* K*exp(-r*tau) - normcdf(-d1)*S_0;
end

function impliedVol = ImpliedVolatilityBlack76Shifted(CP,frwdMarketPrice,K,T,frwdStock,initialVol,shift)
func = @(sigma) (BS_Call_Put_Option_PriceShifted(CP,frwdStock,K,sigma,T,0.0,shift) - frwdMarketPrice).^1.0;
impliedVol = fzero(func,initialVol);

function MakeFigure(X1, YMatrix1, argLegend,titleIn,ylabelTxt)

%CREATEFIGURE(X1,YMATRIX1)
%  X1:  vector of x data
%  YMATRIX1:  matrix of y data

%  Auto-generated by MATLAB on 16-Jan-2012 15:26:40

% Create figure

figure1 = figure('InvertHardcopy','off',...
    'Colormap',[0.061875 0.061875 0.061875;0.06875 0.06875 0.06875;0.075625 0.075625 0.075625;0.0825 0.0825 0.0825;0.089375 0.089375 0.089375;0.09625 0.09625 0.09625;0.103125 0.103125 0.103125;0.11 0.11 0.11;0.146875 0.146875 0.146875;0.18375 0.18375 0.18375;0.220625 0.220625 0.220625;0.2575 0.2575 0.2575;0.294375 0.294375 0.294375;0.33125 0.33125 0.33125;0.368125 0.368125 0.368125;0.405 0.405 0.405;0.441875 0.441875 0.441875;0.47875 0.47875 0.47875;0.515625 0.515625 0.515625;0.5525 0.5525 0.5525;0.589375 0.589375 0.589375;0.62625 0.62625 0.62625;0.663125 0.663125 0.663125;0.7 0.7 0.7;0.711875 0.711875 0.711875;0.72375 0.72375 0.72375;0.735625 0.735625 0.735625;0.7475 0.7475 0.7475;0.759375 0.759375 0.759375;0.77125 0.77125 0.77125;0.783125 0.783125 0.783125;0.795 0.795 0.795;0.806875 0.806875 0.806875;0.81875 0.81875 0.81875;0.830625 0.830625 0.830625;0.8425 0.8425 0.8425;0.854375 0.854375 0.854375;0.86625 0.86625 0.86625;0.878125 0.878125 0.878125;0.89 0.89 0.89;0.853125 0.853125 0.853125;0.81625 0.81625 0.81625;0.779375 0.779375 0.779375;0.7425 0.7425 0.7425;0.705625 0.705625 0.705625;0.66875 0.66875 0.66875;0.631875 0.631875 0.631875;0.595 0.595 0.595;0.558125 0.558125 0.558125;0.52125 0.52125 0.52125;0.484375 0.484375 0.484375;0.4475 0.4475 0.4475;0.410625 0.410625 0.410625;0.37375 0.37375 0.37375;0.336875 0.336875 0.336875;0.3 0.3 0.3;0.28125 0.28125 0.28125;0.2625 0.2625 0.2625;0.24375 0.24375 0.24375;0.225 0.225 0.225;0.20625 0.20625 0.20625;0.1875 0.1875 0.1875;0.16875 0.16875 0.16875;0.15 0.15 0.15],...
    'Color',[1 1 1]);

% Create axes

%axes1 = axes('Parent',figure1,'Color',[1 1 1]);
axes1 = axes('Parent',figure1);
grid on

% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes1,[45 160]);
% Uncomment the following line to preserve the Y-limits of the axes
% ylim(axes1,[19 26]);
% Uncomment the following line to preserve the Z-limits of the axes
% zlim(axes1,[-1 1]);

box(axes1,'on');
hold(axes1,'all');

% Create multiple lines using matrix input to plot
% plot1 = plot(X1,YMatrix1,'Parent',axes1,'MarkerEdgeColor',[0 0 0],...
%     'LineWidth',1,...
%     'Color',[0 0 0]);

plot1 = plot(X1,YMatrix1,'Parent',axes1,...
    'LineWidth',1.5);
set(plot1(1),'Marker','diamond','DisplayName',argLegend{1});
set(plot1(2),'Marker','square','LineStyle','-.',...
    'DisplayName',argLegend{2});
set(plot1(3),'Marker','o','LineStyle','-.','DisplayName',argLegend{3});
set(plot1(4),'DisplayName',argLegend{4});

% Create xlabel

xlabel({'K'});

% Create ylabel

ylabel(ylabelTxt);

% Create title

title(titleIn);

% Create legend

legend1 = legend(axes1,'show');
set(legend1,'Color',[1 1 1]);

